### Election data is taken from: https://data.overheid.nl/community/organization/kiesraad
### And was scraped using the code and procedure taken from here: https://github.com/DIRKMJK/kiesraad

### Groningen

## We begin by processing data for the province of Groningen
## We only annotate the process for Groningen, the process for the remaining provinces is the same

library(openxlsx)
data <- read.csv('combined-gemeente-all-fixed.csv')
eqdata <- read.xlsx('earthquakes_coloured.xlsx')

## Subset election dataset to just Groningen polling stations
groningen <- data[data$contest_name=="Groningen", ]

## Sort the data based on election date
groningen <- groningen[order(groningen$election_date),]

## Remove observations with no location data (LAT is column 11)
groningen <- groningen[complete.cases(groningen[ ,11]),]

## What results after sorting and omitting is
## Obs 1 to 402: 2010 election
## Obs 403 to 821: 2012 election
## Obs 822 to 1224: 2017 election

## The Dutch have their own coordinate system, the Rijksdriehoek, or simply RD system
## The earthquake data uses this system, but I coded the election data in the WGS84 degree decimal system (the more standard system used in GPS and such)
## So some conversion is required
## I have chosen to have our analysis done in RD because RD units are in meters, making distance calculations simpler
## By comparison, WGS84 is in degrees which would require further complicated conversions to become distance

## Function that converts standard WGS84 degree decimal coordinates into the Dutch standard RD system (https://github.com/djvanderlaan/rijksdriehoek/blob/master/R/rijksdriehoek.R)
wgs84_to_rd <- function(lambda, phi) {
  if (!is.numeric(lambda)) stop("lambda needs to be numeric.")
  if (nargs() == 1  && length(lambda) == 2) {
    phi <- lambda[2]
    lambda <- lambda[1]
  } 
  if (!is.numeric(phi)) stop("phi needs to be numeric.")
  if (length(lambda) != length(phi))
    stop("lambda and phi need to have same length.")
  x0 <- 155000.0;
  y0 <- 463000.0;
  phi0 <- 52.15517440;
  lambda0 <-  5.38720621;
  r <- list(
    list(p=0, q=1, r= 190094.945),
    list(p=1, q=1, r= -11832.228),
    list(p=2, q=1, r=   -114.221),
    list(p=0, q=3, r=    -32.391),
    list(p=1, q=0, r=     -0.705),
    list(p=3, q=1, r=     -2.340),
    list(p=1, q=3, r=     -0.608),
    list(p=0, q=2, r=     -0.008),
    list(p=2, q=3, r=      0.148))
  s <- list(
    list(p=1, q=0, s= 309056.544),
    list(p=0, q=2, s=   3638.893),
    list(p=2, q=0, s=     73.077),
    list(p=1, q=2, s=   -157.984),
    list(p=3, q=0, s=     59.788),
    list(p=0, q=1, s=      0.433),
    list(p=2, q=2, s=     -6.439),
    list(p=1, q=1, s=     -0.032),
    list(p=0, q=4, s=      0.092),
    list(p=1, q=4, s=     -0.054))
  dphi = 0.36*(phi - phi0)
  dlambda = 0.36*(lambda - lambda0)
  x <- rep(x0, length(phi))
  y <- rep(y0, length(phi))
  for (i in seq_along(r)) {
    p <- r[[i]][["p"]]
    q <- r[[i]][["q"]]
    ri <- r[[i]][["r"]]
    x <- x +  ri * (dphi^p) * (dlambda^q)
  }
  for (i in seq_along(s)) {
    p <- s[[i]][["p"]]
    q <- s[[i]][["q"]]
    si <- s[[i]][["s"]]
    y <- y +  si * (dphi^p) * (dlambda^q)
  }
  if (nargs() == 1) {
    c(x[1], y[1])
  } else list(x=x, y=y)
}

## Convert standard coordinates to Dutch coordinates, and add as new variables
newcoordg <- wgs84_to_rd(as.numeric(groningen$LON), groningen$LAT)
groningen$RDLON <- newcoordg$x
groningen$RDLAT <- newcoordg$y

## After the 2006 Dutch general elections, the Balkenende gov was installed on 22 Feb 2007
## This means that for the 2010 general election (the oldest election in our dataset), potentially relevant earthquakes are from 22 Feb 2007 to 9 Jun 2010
## This corresponds to earthquakes no. 477 to 659 in the earthquakes dataset

## Correspondingly,
## For the 2012 general elections, the relevant dates from 14 Oct 2010 to 12 Sep 2012, or earthquakes no. 673 to 856
## For the 2017 general elections, the relevant dates from 05 Nov 2012 to 15 Mar 2017, or earthquakes no. 867 to 1368

## Notice that some earthquakes are skipped as they occurred on the dates where negotiations for gov were still in progress

## Coefficients for earthquake PGV calculations, taken from Bommer et al 2017
## We use the PGV_GM column in table 3.1
c1 <- -5.9357
c2 <- 2.4036
c4 <- -1.8819
c4a <- -1.2274
c4b <- -1.7343

## We first create disaggregated earthquake impact variables for earthquakes that occur in different "bins" before the election
## That is, EQs within one year before, EQs more than one year but within two years before, and so on
## Summing up these disaggregated by year values will give us our primary explanatory variable

## This means that for the 2010 general election (the oldest election in our dataset), 
## the first bin of relevant earthquakes are from 10 Jun 2009 to 9 Jun 2010.
## This corresponds to earthquakes no. 620 to 659 in the earthquakes dataset.
## The second bin of relevant earthquakes are from 10 Jun 2008 to 9 Jun 2009.
## This corresponds to earthquakes no. 524 to 619 in the earthquakes dataset.
## The third bin of relevant earthquakes are from 10 Jun 2007 to 9 Jun 2008.
## This corresponds to earthquakes no. 487 to 523 in the earthquakes dataset.
## The fourth bin of relevant earthquakes are from 22 Feb 2007 (date gov installed) to 9 Jun 2007.
## This corresponds to earthquakes no. 477 to 486 in the earthquakes dataset.

## Correspondingly,
## For the 2012 general elections, 
## for the first bin, relevant dates from 11 Sep 2011 to 12 Sep 2012, or earthquakes no. 752 to 856,
## for the second bin, relevant dates from 14 Oct 2010 to 12 Sep 2011 (date gov installed), or earthquakes no. 673 to 751.

## For the 2017 general elections, 
## for the first bin, relevant dates from 14 Mar 2016 to 15 Mar 2017, or earthquakes no. 1246 to 1368
## for the second bin, relevant dates from 14 Mar 2015 to 15 Mar 2015, or earthquakes no. 1126 to 1245
## for the third bin, relevant dates from 14 Mar 2014 to 15 Mar 2015, or earthquakes no. 1038 to 1125
## for the fourth bin, relevant dates from 14 Mar 2013 to 15 Mar 2014, or earthquakes no. 916 to 1037
## for the fifth bin, relevant dates from 05 Nov 2012 to 15 Mar 2013 (date gov installed), or earthquakes no. 867 to 915.

## What follows are several for loops designed to create an impact measure for each observation
## The measure that results from the formulas in Bommer et al is Peak Ground Velocity (PGV), or how fast the ground moves at target location

###########################################################################################

## We will create a main independent variable that simply sums all the PGV in the relevant time periods

groningen$sum_impact_1 <- rep(0)
groningen$sum_impact_2 <- rep(0)
groningen$sum_impact_3 <- rep(0)
groningen$sum_impact_4 <- rep(0)
groningen$sum_impact_5 <- rep(0)

## We've annotated the first loop, the others are the same, just operating on different observations and earthquakes

## 2010 election bin 1
for(j in 1:402){ #outer loop for all observations in the 2010 election
  distance <- rep(NA, 659-619) #create empty vectors to filled by the inner for loop, re-empties it for the next loop
  impact_1 <- rep(NA, 659-619)
  R <- rep(NA, 659-619)
  gR <- rep(NA, 659-619)
  for(i in 620:659){ #inner loop drawing data from the relevant earthquakes
    distance[i-619] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    #pythagoras' theorem to find distance from earthquake epicentre to polling station
    #also converts from meter to km (the earthquake coordinates are in km, the above conversion for polling station coordinates generated meter instead)
    R[i-619] <- sqrt(distance[i-619]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2) #equation 3.2 from the Bommer report
    if(R[i-619] < 6.32){ #an ifelse statement that applies the appropriate equation 3.3 depending on the calculated distance
      gR[i-619] <- c4 * log(R[i-619])
    } else if(R[i-619] < 11.62){
      gR[i-619] <- c4 * log(6.32) + c4a * log(R[i-619]/6.32)
    } else{
      gR[i-619] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-619]/11.62)
    }
    impact_1[i-619] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-619]) #equation 3.1, generates the PGV
  }
  groningen$sum_impact_1[j] <- sum(impact_1)  #sum of all earthquake impact in the relevant time period for that polling station
}

## 2010 election bin 2
for(j in 1:402){
  distance <- rep(NA, 619-523)
  impact_2 <- rep(NA, 619-523)
  R <- rep(NA, 619-523)
  gR <- rep(NA, 619-523)
  for(i in 524:619){
    distance[i-523] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-523] <- sqrt(distance[i-523]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-523] < 6.32){
      gR[i-523] <- c4 * log(R[i-523])
    } else if(R[i-523] < 11.62){
      gR[i-523] <- c4 * log(6.32) + c4a * log(R[i-523]/6.32)
    } else{
      gR[i-523] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-523]/11.62)
    }
    impact_2[i-523] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-523])
  }
  groningen$sum_impact_2[j] <- sum(impact_2)
}

## 2010 election bin 3
for(j in 1:402){
  distance <- rep(NA, 523-486)
  impact_3 <- rep(NA, 523-486)
  R <- rep(NA, 523-486)
  gR <- rep(NA, 523-486)
  for(i in 487:523){
    distance[i-486] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-486] <- sqrt(distance[i-486]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-486] < 6.32){
      gR[i-486] <- c4 * log(R[i-486])
    } else if(R[i-486] < 11.62){
      gR[i-486] <- c4 * log(6.32) + c4a * log(R[i-486]/6.32)
    } else{
      gR[i-486] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-486]/11.62)
    }
    impact_3[i-486] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-486])
  }
  groningen$sum_impact_3[j] <- sum(impact_3)
}

## 2010 election bin 4
for(j in 1:402){
  distance <- rep(NA, 486-476)
  impact_4 <- rep(NA, 486-476)
  R <- rep(NA, 486-476)
  gR <- rep(NA, 486-476)
  for(i in 477:486){
    distance[i-476] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-476] <- sqrt(distance[i-476]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-476] < 6.32){
      gR[i-476] <- c4 * log(R[i-476])
    } else if(R[i-476] < 11.62){
      gR[i-476] <- c4 * log(6.32) + c4a * log(R[i-476]/6.32)
    } else{
      gR[i-476] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-476]/11.62)
    }
    impact_4[i-476] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-476])
  }
  groningen$sum_impact_4[j] <- sum(impact_4)
}

## 2012 election bin 1
for(j in 403:821){
  distance <- rep(NA, 856-751)
  impact_1 <- rep(NA, 856-751)
  R <- rep(NA, 856-751)
  gR <- rep(NA, 856-751)
  for(i in 752:856){
    distance[i-751] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-751] <- sqrt(distance[i-751]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-751] < 6.32){
      gR[i-751] <- c4 * log(R[i-751])
    } else if(R[i-751] < 11.62){
      gR[i-751] <- c4 * log(6.32) + c4a * log(R[i-751]/6.32)
    } else{
      gR[i-751] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-751]/11.62)
    }
    impact_1[i-751] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-751])
  }
  groningen$sum_impact_1[j] <- sum(impact_1)
}

## 2012 election bin 2
for(j in 403:821){
  distance <- rep(NA, 751-673)
  impact_2 <- rep(NA, 751-673)
  R <- rep(NA, 751-673)
  gR <- rep(NA, 751-673)
  for(i in 674:751){
    distance[i-673] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-673] <- sqrt(distance[i-673]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-673] < 6.32){
      gR[i-673] <- c4 * log(R[i-673])
    } else if(R[i-673] < 11.62){
      gR[i-673] <- c4 * log(6.32) + c4a * log(R[i-673]/6.32)
    } else{
      gR[i-673] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-673]/11.62)
    }
    impact_2[i-673] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-673])
  }
  groningen$sum_impact_2[j] <- sum(impact_2)
}

## 2017 election bin 1
for(j in 822:1224){
  distance <- rep(NA, 1368-1245)
  impact_1 <- rep(NA, 1368-1245)
  R <- rep(NA, 1368-1245)
  gR <- rep(NA, 1368-1245)
  for(i in 1246:1368){
    distance[i-1245] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-1245] <- sqrt(distance[i-1245]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1245] < 6.32){
      gR[i-1245] <- c4 * log(R[i-1245])
    } else if(R[i-1245] < 11.62){
      gR[i-1245] <- c4 * log(6.32) + c4a * log(R[i-1245]/6.32)
    } else{
      gR[i-1245] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1245]/11.62)
    }
    impact_1[i-1245] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1245])
  }
  groningen$sum_impact_1[j] <- sum(impact_1)
}

## 2017 election bin 2
for(j in 822:1224){
  distance <- rep(NA, 1245-1125)
  impact_2 <- rep(NA, 1245-1125)
  R <- rep(NA, 1245-1125)
  gR <- rep(NA, 1245-1125)
  for(i in 1126:1245){
    distance[i-1125] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-1125] <- sqrt(distance[i-1125]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1125] < 6.32){
      gR[i-1125] <- c4 * log(R[i-1125])
    } else if(R[i-1125] < 11.62){
      gR[i-1125] <- c4 * log(6.32) + c4a * log(R[i-1125]/6.32)
    } else{
      gR[i-1125] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1125]/11.62)
    }
    impact_2[i-1125] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1125])
  }
  groningen$sum_impact_2[j] <- sum(impact_2)
}

## 2017 election bin 3
for(j in 822:1224){
  distance <- rep(NA, 1125-1037)
  impact_3 <- rep(NA, 1125-1037)
  R <- rep(NA, 1125-1037)
  gR <- rep(NA, 1125-1037)
  for(i in 1038:1125){
    distance[i-1037] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-1037] <- sqrt(distance[i-1037]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1037] < 6.32){
      gR[i-1037] <- c4 * log(R[i-1037])
    } else if(R[i-1037] < 11.62){
      gR[i-1037] <- c4 * log(6.32) + c4a * log(R[i-1037]/6.32)
    } else{
      gR[i-1037] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1037]/11.62)
    }
    impact_3[i-1037] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1037])
  }
  groningen$sum_impact_3[j] <- sum(impact_3)
}

## 2017 election bin 4
for(j in 822:1224){
  distance <- rep(NA, 1037-915)
  impact_4 <- rep(NA, 1037-915)
  R <- rep(NA, 1037-915)
  gR <- rep(NA, 1037-915)
  for(i in 916:1037){
    distance[i-915] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-915] <- sqrt(distance[i-915]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-915] < 6.32){
      gR[i-915] <- c4 * log(R[i-915])
    } else if(R[i-915] < 11.62){
      gR[i-915] <- c4 * log(6.32) + c4a * log(R[i-915]/6.32)
    } else{
      gR[i-915] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-915]/11.62)
    }
    impact_4[i-915] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-915])
  }
  groningen$sum_impact_4[j] <- sum(impact_4)
}

## 2017 election bin 5
for(j in 822:1224){
  distance <- rep(NA, 915-866)
  impact_5 <- rep(NA, 915-866)
  R <- rep(NA, 915-866)
  gR <- rep(NA, 915-866)
  for(i in 867:915){
    distance[i-866] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-866] <- sqrt(distance[i-866]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-866] < 6.32){
      gR[i-866] <- c4 * log(R[i-866])
    } else if(R[i-866] < 11.62){
      gR[i-866] <- c4 * log(6.32) + c4a * log(R[i-866]/6.32)
    } else{
      gR[i-866] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-866]/11.62)
    }
    impact_5[i-866] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-866])
  }
  groningen$sum_impact_5[j] <- sum(impact_5)
}

## Summing the 5 bins

groningen$sum_impact <- groningen$sum_impact_1 + groningen$sum_impact_2 + groningen$sum_impact_3 + groningen$sum_impact_4 + groningen$sum_impact_5

###########################################################################################

## We repeat all the above processes for the provinces of Drenthe and Friesland

### Drenthe

drenthe <- data[data$contest_name=="Assen", ]

drenthe <- drenthe[order(drenthe$election_date),]

drenthe <- drenthe[complete.cases(drenthe[ ,11]),]

## What results after sorting and omitting is
## Obs 1 to 334: 2010 election
## Obs 335 to 671: 2012 election
## Obs 672 to 1012: 2017 election

newcoordd <- wgs84_to_rd(as.numeric(drenthe$LON), drenthe$LAT)
drenthe$RDLON <- newcoordd$x
drenthe$RDLAT <- newcoordd$y

drenthe$sum_impact_1 <- rep(0)
drenthe$sum_impact_2 <- rep(0)
drenthe$sum_impact_3 <- rep(0)
drenthe$sum_impact_4 <- rep(0)
drenthe$sum_impact_5 <- rep(0)

## 2010 election bin 1
for(j in 1:334){
  distance <- rep(NA, 659-619)
  impact_1 <- rep(NA, 659-619)
  R <- rep(NA, 659-619)
  gR <- rep(NA, 659-619)
  for(i in 620:659){
    distance[i-619] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-619] <- sqrt(distance[i-619]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-619] < 6.32){
      gR[i-619] <- c4 * log(R[i-619])
    } else if(R[i-619] < 11.62){
      gR[i-619] <- c4 * log(6.32) + c4a * log(R[i-619]/6.32)
    } else{
      gR[i-619] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-619]/11.62)
    }
    impact_1[i-619] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-619])
  }
  drenthe$sum_impact_1[j] <- sum(impact_1)
}

## 2010 election bin 2
for(j in 1:334){
  distance <- rep(NA, 619-523)
  impact_2 <- rep(NA, 619-523)
  R <- rep(NA, 619-523)
  gR <- rep(NA, 619-523)
  for(i in 524:619){
    distance[i-523] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-523] <- sqrt(distance[i-523]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-523] < 6.32){
      gR[i-523] <- c4 * log(R[i-523])
    } else if(R[i-523] < 11.62){
      gR[i-523] <- c4 * log(6.32) + c4a * log(R[i-523]/6.32)
    } else{
      gR[i-523] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-523]/11.62)
    }
    impact_2[i-523] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-523])
  }
  drenthe$sum_impact_2[j] <- sum(impact_2)
}

## 2010 election bin 3
for(j in 1:334){
  distance <- rep(NA, 523-486)
  impact_3 <- rep(NA, 523-486)
  R <- rep(NA, 523-486)
  gR <- rep(NA, 523-486)
  for(i in 487:523){
    distance[i-486] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-486] <- sqrt(distance[i-486]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-486] < 6.32){
      gR[i-486] <- c4 * log(R[i-486])
    } else if(R[i-486] < 11.62){
      gR[i-486] <- c4 * log(6.32) + c4a * log(R[i-486]/6.32)
    } else{
      gR[i-486] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-486]/11.62)
    }
    impact_3[i-486] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-486])
  }
  drenthe$sum_impact_3[j] <- sum(impact_3)
}

## 2010 election bin 4
for(j in 1:334){
  distance <- rep(NA, 486-476)
  impact_4 <- rep(NA, 486-476)
  R <- rep(NA, 486-476)
  gR <- rep(NA, 486-476)
  for(i in 477:486){
    distance[i-476] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-476] <- sqrt(distance[i-476]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-476] < 6.32){
      gR[i-476] <- c4 * log(R[i-476])
    } else if(R[i-476] < 11.62){
      gR[i-476] <- c4 * log(6.32) + c4a * log(R[i-476]/6.32)
    } else{
      gR[i-476] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-476]/11.62)
    }
    impact_4[i-476] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-476])
  }
  drenthe$sum_impact_4[j] <- sum(impact_4)
}

## 2012 election bin 1
for(j in 335:671){
  distance <- rep(NA, 856-751)
  impact_1 <- rep(NA, 856-751)
  R <- rep(NA, 856-751)
  gR <- rep(NA, 856-751)
  for(i in 752:856){
    distance[i-751] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-751] <- sqrt(distance[i-751]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-751] < 6.32){
      gR[i-751] <- c4 * log(R[i-751])
    } else if(R[i-751] < 11.62){
      gR[i-751] <- c4 * log(6.32) + c4a * log(R[i-751]/6.32)
    } else{
      gR[i-751] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-751]/11.62)
    }
    impact_1[i-751] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-751])
  }
  drenthe$sum_impact_1[j] <- sum(impact_1)
}

## 2012 election bin 2
for(j in 335:671){
  distance <- rep(NA, 751-673)
  impact_2 <- rep(NA, 751-673)
  R <- rep(NA, 751-673)
  gR <- rep(NA, 751-673)
  for(i in 674:751){
    distance[i-673] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-673] <- sqrt(distance[i-673]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-673] < 6.32){
      gR[i-673] <- c4 * log(R[i-673])
    } else if(R[i-673] < 11.62){
      gR[i-673] <- c4 * log(6.32) + c4a * log(R[i-673]/6.32)
    } else{
      gR[i-673] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-673]/11.62)
    }
    impact_2[i-673] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-673])
  }
  drenthe$sum_impact_2[j] <- sum(impact_2)
}

## 2017 election bin 1
for(j in 672:1012){
  distance <- rep(NA, 1368-1245)
  impact_1 <- rep(NA, 1368-1245)
  R <- rep(NA, 1368-1245)
  gR <- rep(NA, 1368-1245)
  for(i in 1246:1368){
    distance[i-1245] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-1245] <- sqrt(distance[i-1245]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1245] < 6.32){
      gR[i-1245] <- c4 * log(R[i-1245])
    } else if(R[i-1245] < 11.62){
      gR[i-1245] <- c4 * log(6.32) + c4a * log(R[i-1245]/6.32)
    } else{
      gR[i-1245] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1245]/11.62)
    }
    impact_1[i-1245] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1245])
  }
  drenthe$sum_impact_1[j] <- sum(impact_1)
}

## 2017 election bin 2
for(j in 672:1012){
  distance <- rep(NA, 1245-1125)
  impact_2 <- rep(NA, 1245-1125)
  R <- rep(NA, 1245-1125)
  gR <- rep(NA, 1245-1125)
  for(i in 1126:1245){
    distance[i-1125] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-1125] <- sqrt(distance[i-1125]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1125] < 6.32){
      gR[i-1125] <- c4 * log(R[i-1125])
    } else if(R[i-1125] < 11.62){
      gR[i-1125] <- c4 * log(6.32) + c4a * log(R[i-1125]/6.32)
    } else{
      gR[i-1125] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1125]/11.62)
    }
    impact_2[i-1125] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1125])
  }
  drenthe$sum_impact_2[j] <- sum(impact_2)
}

## 2017 election bin 3
for(j in 672:1012){
  distance <- rep(NA, 1125-1037)
  impact_3 <- rep(NA, 1125-1037)
  R <- rep(NA, 1125-1037)
  gR <- rep(NA, 1125-1037)
  for(i in 1038:1125){
    distance[i-1037] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-1037] <- sqrt(distance[i-1037]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1037] < 6.32){
      gR[i-1037] <- c4 * log(R[i-1037])
    } else if(R[i-1037] < 11.62){
      gR[i-1037] <- c4 * log(6.32) + c4a * log(R[i-1037]/6.32)
    } else{
      gR[i-1037] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1037]/11.62)
    }
    impact_3[i-1037] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1037])
  }
  drenthe$sum_impact_3[j] <- sum(impact_3)
}

## 2017 election bin 4
for(j in 672:1012){
  distance <- rep(NA, 1037-915)
  impact_4 <- rep(NA, 1037-915)
  R <- rep(NA, 1037-915)
  gR <- rep(NA, 1037-915)
  for(i in 916:1037){
    distance[i-915] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-915] <- sqrt(distance[i-915]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-915] < 6.32){
      gR[i-915] <- c4 * log(R[i-915])
    } else if(R[i-915] < 11.62){
      gR[i-915] <- c4 * log(6.32) + c4a * log(R[i-915]/6.32)
    } else{
      gR[i-915] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-915]/11.62)
    }
    impact_4[i-915] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-915])
  }
  drenthe$sum_impact_4[j] <- sum(impact_4)
}

## 2017 election bin 5
for(j in 672:1012){
  distance <- rep(NA, 915-866)
  impact_5 <- rep(NA, 915-866)
  R <- rep(NA, 915-866)
  gR <- rep(NA, 915-866)
  for(i in 867:915){
    distance[i-866] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-866] <- sqrt(distance[i-866]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-866] < 6.32){
      gR[i-866] <- c4 * log(R[i-866])
    } else if(R[i-866] < 11.62){
      gR[i-866] <- c4 * log(6.32) + c4a * log(R[i-866]/6.32)
    } else{
      gR[i-866] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-866]/11.62)
    }
    impact_5[i-866] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-866])
  }
  drenthe$sum_impact_5[j] <- sum(impact_5)
}

## Summing the 5 bins

drenthe$sum_impact <- drenthe$sum_impact_1 + drenthe$sum_impact_2 + drenthe$sum_impact_3 + drenthe$sum_impact_4 + drenthe$sum_impact_5

###########################################################################################

### Friesland

friesland <- data[data$contest_name=="Leeuwarden", ]

friesland <- friesland[order(friesland$election_date),]

friesland <- friesland[complete.cases(friesland[ ,11]),]

## What results after sorting and omitting is
## Obs 1 to 436: 2010 election
## Obs 437 to 895: 2012 election
## Obs 896 to 1337: 2017 election

newcoordf <- wgs84_to_rd(as.numeric(friesland$LON), friesland$LAT)
friesland$RDLON <- newcoordf$x
friesland$RDLAT <- newcoordf$y

friesland$sum_impact_1 <- rep(0)
friesland$sum_impact_2 <- rep(0)
friesland$sum_impact_3 <- rep(0)
friesland$sum_impact_4 <- rep(0)
friesland$sum_impact_5 <- rep(0)

## 2010 election bin 1
for(j in 1:436){
  distance <- rep(NA, 659-619)
  impact_1 <- rep(NA, 659-619)
  R <- rep(NA, 659-619)
  gR <- rep(NA, 659-619)
  for(i in 620:659){
    distance[i-619] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-619] <- sqrt(distance[i-619]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-619] < 6.32){
      gR[i-619] <- c4 * log(R[i-619])
    } else if(R[i-619] < 11.62){
      gR[i-619] <- c4 * log(6.32) + c4a * log(R[i-619]/6.32)
    } else{
      gR[i-619] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-619]/11.62)
    }
    impact_1[i-619] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-619])
  }
  friesland$sum_impact_1[j] <- sum(impact_1)
}

## 2010 election bin 2
for(j in 1:436){
  distance <- rep(NA, 619-523)
  impact_2 <- rep(NA, 619-523)
  R <- rep(NA, 619-523)
  gR <- rep(NA, 619-523)
  for(i in 524:619){
    distance[i-523] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-523] <- sqrt(distance[i-523]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-523] < 6.32){
      gR[i-523] <- c4 * log(R[i-523])
    } else if(R[i-523] < 11.62){
      gR[i-523] <- c4 * log(6.32) + c4a * log(R[i-523]/6.32)
    } else{
      gR[i-523] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-523]/11.62)
    }
    impact_2[i-523] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-523])
  }
  friesland$sum_impact_2[j] <- sum(impact_2)
}

## 2010 election bin 3
for(j in 1:436){
  distance <- rep(NA, 523-486)
  impact_3 <- rep(NA, 523-486)
  R <- rep(NA, 523-486)
  gR <- rep(NA, 523-486)
  for(i in 487:523){
    distance[i-486] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-486] <- sqrt(distance[i-486]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-486] < 6.32){
      gR[i-486] <- c4 * log(R[i-486])
    } else if(R[i-486] < 11.62){
      gR[i-486] <- c4 * log(6.32) + c4a * log(R[i-486]/6.32)
    } else{
      gR[i-486] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-486]/11.62)
    }
    impact_3[i-486] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-486])
  }
  friesland$sum_impact_3[j] <- sum(impact_3)
}

## 2010 election bin 4
for(j in 1:436){
  distance <- rep(NA, 486-476)
  impact_4 <- rep(NA, 486-476)
  R <- rep(NA, 486-476)
  gR <- rep(NA, 486-476)
  for(i in 477:486){
    distance[i-476] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-476] <- sqrt(distance[i-476]^2 + (exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083))^2)
    if(R[i-476] < 6.32){
      gR[i-476] <- c4 * log(R[i-476])
    } else if(R[i-476] < 11.62){
      gR[i-476] <- c4 * log(6.32) + c4a * log(R[i-476]/6.32)
    } else{
      gR[i-476] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-476]/11.62)
    }
    impact_4[i-476] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-476])
  }
  friesland$sum_impact_4[j] <- sum(impact_4)
}

## 2012 election bin 1
for(j in 437:895){
  distance <- rep(NA, 856-751)
  impact_1 <- rep(NA, 856-751)
  R <- rep(NA, 856-751)
  gR <- rep(NA, 856-751)
  for(i in 752:856){
    distance[i-751] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-751] <- sqrt(distance[i-751]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-751] < 6.32){
      gR[i-751] <- c4 * log(R[i-751])
    } else if(R[i-751] < 11.62){
      gR[i-751] <- c4 * log(6.32) + c4a * log(R[i-751]/6.32)
    } else{
      gR[i-751] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-751]/11.62)
    }
    impact_1[i-751] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-751])
  }
  friesland$sum_impact_1[j] <- sum(impact_1)
}

## 2012 election bin 2
for(j in 437:895){
  distance <- rep(NA, 751-673)
  impact_2 <- rep(NA, 751-673)
  R <- rep(NA, 751-673)
  gR <- rep(NA, 751-673)
  for(i in 674:751){
    distance[i-673] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-673] <- sqrt(distance[i-673]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-673] < 6.32){
      gR[i-673] <- c4 * log(R[i-673])
    } else if(R[i-673] < 11.62){
      gR[i-673] <- c4 * log(6.32) + c4a * log(R[i-673]/6.32)
    } else{
      gR[i-673] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-673]/11.62)
    }
    impact_2[i-673] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-673])
  }
  friesland$sum_impact_2[j] <- sum(impact_2)
}

## 2017 election bin 1
for(j in 896:1337){
  distance <- rep(NA, 1368-1245)
  impact_1 <- rep(NA, 1368-1245)
  R <- rep(NA, 1368-1245)
  gR <- rep(NA, 1368-1245)
  for(i in 1246:1368){
    distance[i-1245] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-1245] <- sqrt(distance[i-1245]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1245] < 6.32){
      gR[i-1245] <- c4 * log(R[i-1245])
    } else if(R[i-1245] < 11.62){
      gR[i-1245] <- c4 * log(6.32) + c4a * log(R[i-1245]/6.32)
    } else{
      gR[i-1245] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1245]/11.62)
    }
    impact_1[i-1245] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1245])
  }
  friesland$sum_impact_1[j] <- sum(impact_1)
}

## 2017 election bin 2
for(j in 896:1337){
  distance <- rep(NA, 1245-1125)
  impact_2 <- rep(NA, 1245-1125)
  R <- rep(NA, 1245-1125)
  gR <- rep(NA, 1245-1125)
  for(i in 1126:1245){
    distance[i-1125] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-1125] <- sqrt(distance[i-1125]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1125] < 6.32){
      gR[i-1125] <- c4 * log(R[i-1125])
    } else if(R[i-1125] < 11.62){
      gR[i-1125] <- c4 * log(6.32) + c4a * log(R[i-1125]/6.32)
    } else{
      gR[i-1125] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1125]/11.62)
    }
    impact_2[i-1125] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1125])
  }
  friesland$sum_impact_2[j] <- sum(impact_2)
}

## 2017 election bin 3
for(j in 896:1337){
  distance <- rep(NA, 1125-1037)
  impact_3 <- rep(NA, 1125-1037)
  R <- rep(NA, 1125-1037)
  gR <- rep(NA, 1125-1037)
  for(i in 1038:1125){
    distance[i-1037] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-1037] <- sqrt(distance[i-1037]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-1037] < 6.32){
      gR[i-1037] <- c4 * log(R[i-1037])
    } else if(R[i-1037] < 11.62){
      gR[i-1037] <- c4 * log(6.32) + c4a * log(R[i-1037]/6.32)
    } else{
      gR[i-1037] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-1037]/11.62)
    }
    impact_3[i-1037] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-1037])
  }
  friesland$sum_impact_3[j] <- sum(impact_3)
}

## 2017 election bin 4
for(j in 896:1337){
  distance <- rep(NA, 1037-915)
  impact_4 <- rep(NA, 1037-915)
  R <- rep(NA, 1037-915)
  gR <- rep(NA, 1037-915)
  for(i in 916:1037){
    distance[i-915] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-915] <- sqrt(distance[i-915]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-915] < 6.32){
      gR[i-915] <- c4 * log(R[i-915])
    } else if(R[i-915] < 11.62){
      gR[i-915] <- c4 * log(6.32) + c4a * log(R[i-915]/6.32)
    } else{
      gR[i-915] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-915]/11.62)
    }
    impact_4[i-915] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-915])
  }
  friesland$sum_impact_4[j] <- sum(impact_4)
}

## 2017 election bin 5
for(j in 896:1337){
  distance <- rep(NA, 915-866)
  impact_5 <- rep(NA, 915-866)
  R <- rep(NA, 915-866)
  gR <- rep(NA, 915-866)
  for(i in 867:915){
    distance[i-866] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-866] <- sqrt(distance[i-866]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-866] < 6.32){
      gR[i-866] <- c4 * log(R[i-866])
    } else if(R[i-866] < 11.62){
      gR[i-866] <- c4 * log(6.32) + c4a * log(R[i-866]/6.32)
    } else{
      gR[i-866] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-866]/11.62)
    }
    impact_5[i-866] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-866])
  }
  friesland$sum_impact_5[j] <- sum(impact_5)
}

## Summing the 5 bins

friesland$sum_impact <- friesland$sum_impact_1 + friesland$sum_impact_2 + friesland$sum_impact_3 + friesland$sum_impact_4 + friesland$sum_impact_5

###########################################################################################

### Here we create data for our "placebo test" models using prior elections

## Groningen

groningen$sum_impact_p <- rep(NA)

## What results after sorting and omitting is
## Obs 1 to 402: 2010 election
## Obs 403 to 821: 2012 election

## For the 2012 general elections, the relevant dates from 14 Oct 2010 to 12 Sep 2012, or earthquakes no. 673 to 856
## For the 2017 general elections, the relevant dates from 05 Nov 2012 to 15 Mar 2017, or earthquakes no. 867 to 1368

## 2010 election
for(j in 1:402){ 
  distance <- rep(NA, 856-672)
  impact_p <- rep(NA, 856-672)
  R <- rep(NA, 856-672)
  gR <- rep(NA, 856-672)
  for(i in 673:856){
    distance[i-672] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-672] <- sqrt(distance[i-672]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-672] < 6.32){
      gR[i-672] <- c4 * log(R[i-672])
    } else if(R[i-672] < 11.62){
      gR[i-672] <- c4 * log(6.32) + c4a * log(R[i-672]/6.32)
    } else{
      gR[i-672] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-672]/11.62)
    }
    impact_p[i-672] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-672])
  }
  groningen$sum_impact_p[j] <- sum(impact_p)
}

## 2012 election
for(j in 403:821){
  distance <- rep(NA, 1368-866)
  impact_p <- rep(NA, 1368-866)
  R <- rep(NA, 1368-866)
  gR <- rep(NA, 1368-866)
  for(i in 867:1368){
    distance[i-866] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-groningen$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-groningen$RDLON[j])^2) / 1000
    R[i-866] <- sqrt(distance[i-866]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-866] < 6.32){
      gR[i-866] <- c4 * log(R[i-866])
    } else if(R[i-866] < 11.62){
      gR[i-866] <- c4 * log(6.32) + c4a * log(R[i-866]/6.32)
    } else{
      gR[i-866] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-866]/11.62)
    }
    impact_p[i-866] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-866])
  }
  groningen$sum_impact_p[j] <- sum(impact_p)
}

###########################################################################################

## Drenthe

drenthe$sum_impact_p <- rep(NA)

## What results after sorting and omitting is
## Obs 1 to 334: 2010 election
## Obs 335 to 671: 2012 election

## For the 2012 general elections, the relevant dates from 14 Oct 2010 to 12 Sep 2012, or earthquakes no. 673 to 856
## For the 2017 general elections, the relevant dates from 05 Nov 2012 to 15 Mar 2017, or earthquakes no. 867 to 1368

## 2010 election
for(j in 1:334){ 
  distance <- rep(NA, 856-672)
  impact_p <- rep(NA, 856-672)
  R <- rep(NA, 856-672)
  gR <- rep(NA, 856-672)
  for(i in 673:856){
    distance[i-672] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-672] <- sqrt(distance[i-672]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-672] < 6.32){
      gR[i-672] <- c4 * log(R[i-672])
    } else if(R[i-672] < 11.62){
      gR[i-672] <- c4 * log(6.32) + c4a * log(R[i-672]/6.32)
    } else{
      gR[i-672] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-672]/11.62)
    }
    impact_p[i-672] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-672])
  }
  drenthe$sum_impact_p[j] <- sum(impact_p)
}

## 2012 election
for(j in 335:671){
  distance <- rep(NA, 1368-866)
  impact_p <- rep(NA, 1368-866)
  R <- rep(NA, 1368-866)
  gR <- rep(NA, 1368-866)
  for(i in 867:1368){
    distance[i-866] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-drenthe$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-drenthe$RDLON[j])^2) / 1000
    R[i-866] <- sqrt(distance[i-866]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-866] < 6.32){
      gR[i-866] <- c4 * log(R[i-866])
    } else if(R[i-866] < 11.62){
      gR[i-866] <- c4 * log(6.32) + c4a * log(R[i-866]/6.32)
    } else{
      gR[i-866] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-866]/11.62)
    }
    impact_p[i-866] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-866])
  }
  drenthe$sum_impact_p[j] <- sum(impact_p)
}

###########################################################################################

## Friesland

friesland$sum_impact_p <- rep(NA)

## What results after sorting and omitting is
## Obs 1 to 436: 2010 election
## Obs 437 to 895: 2012 election

## For the 2012 general elections, the relevant dates from 14 Oct 2010 to 12 Sep 2012, or earthquakes no. 673 to 856
## For the 2017 general elections, the relevant dates from 05 Nov 2012 to 15 Mar 2017, or earthquakes no. 867 to 1368

## 2010 election
for(j in 1:436){ 
  distance <- rep(NA, 856-672)
  impact_p <- rep(NA, 856-672)
  R <- rep(NA, 856-672)
  gR <- rep(NA, 856-672)
  for(i in 673:856){
    distance[i-672] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-672] <- sqrt(distance[i-672]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-672] < 6.32){
      gR[i-672] <- c4 * log(R[i-672])
    } else if(R[i-672] < 11.62){
      gR[i-672] <- c4 * log(6.32) + c4a * log(R[i-672]/6.32)
    } else{
      gR[i-672] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-672]/11.62)
    }
    impact_p[i-672] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-672])
  }
  friesland$sum_impact_p[j] <- sum(impact_p)
}

## 2012 election
for(j in 437:895){
  distance <- rep(NA, 1368-866)
  impact_p <- rep(NA, 1368-866)
  R <- rep(NA, 1368-866)
  gR <- rep(NA, 1368-866)
  for(i in 867:1368){
    distance[i-866] <- sqrt((eqdata$LAT[eqdata$EQ_ID==i]-friesland$RDLAT[j])^2 + (eqdata$LON[eqdata$EQ_ID==i]-friesland$RDLON[j])^2) / 1000
    R[i-866] <- sqrt(distance[i-866]^2 + exp(0.4233*eqdata$MAG[eqdata$EQ_ID==i] - 0.6083)^2)
    if(R[i-866] < 6.32){
      gR[i-866] <- c4 * log(R[i-866])
    } else if(R[i-866] < 11.62){
      gR[i-866] <- c4 * log(6.32) + c4a * log(R[i-866]/6.32)
    } else{
      gR[i-866] <- c4 * log(6.32) + c4a * log(11.62/6.32) + c4b * log(R[i-866]/11.62)
    }
    impact_p[i-866] <- exp(c1 + c2*eqdata$MAG[eqdata$EQ_ID==i] + gR[i-866])
  }
  friesland$sum_impact_p[j] <- sum(impact_p)
}

  ### Remove 2021 election data (don't have earthquake data to match)
  
  ## Groningen
  groningen <- groningen[groningen$election_id != "TK2021",]
  
  ## Drenthe
  drenthe <- drenthe[drenthe$election_id != "TK2021",]
  
  ## Friesland
  friesland <- friesland[friesland$election_id != "TK2021",]
  
  ### We now combine the data from the three provinces
  
  ## All
  
  allprovinces <- rbind(groningen, drenthe, friesland)
  
  ## Create a variable that is vote share for GroenLinks, the major Green party
  allprovinces$GL_pct <- allprovinces$GROENLINKS / allprovinces$total_counted
  
  ## Drill distance control
  
  #RD coordinates of active gas wells
  #t Zandt	248027.62	600458.4
  #Amsweer	256193	591550
  #Bierum	254705	599497
  #de Eeker	259676	577292
  #de Paauwen	246111.77	588405.41
  #Eemskanaal	241602.09	584344.27
  #Froombosch	248302	578964
  #Kooipolder	246879	580965
  #Leermens	249897.88	596901.72
  #Oudeweg	256116.55	585601.04
  #Overschild	250592	590664
  #Sappemeer	249591	575459
  #Schaapbulten	257607.95	588337.01
  #Scheemderzwag	257064	578126
  #Siddeburen	253054	587461
  #Slochteren	246347	579286
  #Spitsbergen	252380	577164
  #Ten Post	245707	591156
  #Tjuchem	254940.35	588055.22
  #Tusschenklappen	254217.54	575229.01
  #Zuiderpolder	261973.21	580692.88
  #Zuiderveen	264336	573311
  
  allprovinces$drill <- (sqrt((allprovinces$RDLON - 248027.62)^2 + (allprovinces$RDLAT - 600458.4)^2) + 
                           sqrt((allprovinces$RDLON - 256193)^2 + (allprovinces$RDLAT - 591550)^2) + 
                           sqrt((allprovinces$RDLON - 254705)^2 + (allprovinces$RDLAT - 599497)^2) + 
                           sqrt((allprovinces$RDLON - 259676)^2 + (allprovinces$RDLAT - 577292)^2) + 
                           sqrt((allprovinces$RDLON - 246111.77)^2 + (allprovinces$RDLAT - 588405.41)^2) + 
                           sqrt((allprovinces$RDLON - 241602.09)^2 + (allprovinces$RDLAT - 584344.27)^2) + 
                           sqrt((allprovinces$RDLON - 248302)^2 + (allprovinces$RDLAT - 578964)^2) + 
                           sqrt((allprovinces$RDLON - 246879)^2 + (allprovinces$RDLAT - 580965)^2) + 
                           sqrt((allprovinces$RDLON - 249897.88)^2 + (allprovinces$RDLAT - 596901.72)^2) + 
                           sqrt((allprovinces$RDLON - 256116.55)^2 + (allprovinces$RDLAT - 585601.04)^2) + 
                           sqrt((allprovinces$RDLON - 250592)^2 + (allprovinces$RDLAT - 590664)^2) + 
                           sqrt((allprovinces$RDLON - 249591)^2 + (allprovinces$RDLAT - 575459)^2) + 
                           sqrt((allprovinces$RDLON - 257607.95)^2 + (allprovinces$RDLAT - 588337.01)^2) + 
                           sqrt((allprovinces$RDLON - 257064)^2 + (allprovinces$RDLAT - 578126)^2) + 
                           sqrt((allprovinces$RDLON - 253054)^2 + (allprovinces$RDLAT - 587461)^2) + 
                           sqrt((allprovinces$RDLON - 246347)^2 + (allprovinces$RDLAT - 579286)^2) + 
                           sqrt((allprovinces$RDLON - 252380)^2 + (allprovinces$RDLAT - 577164)^2) + 
                           sqrt((allprovinces$RDLON - 245707)^2 + (allprovinces$RDLAT - 591156)^2) + 
                           sqrt((allprovinces$RDLON - 254940.35)^2 + (allprovinces$RDLAT - 588055.22)^2) + 
                           sqrt((allprovinces$RDLON - 254217.54)^2 + (allprovinces$RDLAT - 575229.01)^2) + 
                           sqrt((allprovinces$RDLON - 261973.21)^2 + (allprovinces$RDLAT - 580692.88)^2) + 
                           sqrt((allprovinces$RDLON - 264336)^2 + (allprovinces$RDLAT - 573311)^2))/22
  
  ## Subset the data, for 4 bins regression
  
  bin4 <- subset(allprovinces, election_id != "TK2012")
  
  ### Run overall regressions
  
  ## Regressions using all provinces, main earthquake variable
  
  library(sandwich)
  library(lmtest)
  all_sumCD <- lm(GL_pct ~ sum_impact + cast + drill, allprovinces)
  all_sumCD_clustered <- coeftest(all_sumCD, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_clusterSE <- all_sumCD_clustered[,2]
  all_sumCDE <- lm(GL_pct ~ sum_impact + cast + drill + as.factor(election_id), allprovinces)
  all_sumCDE_clustered <- coeftest(all_sumCDE, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCDE_clusterSE <- all_sumCDE_clustered[,2]
  all_sumACD <- lm(GL_pct ~ sum_impact + cast + drill + as.factor(managing_authority), allprovinces)
  all_sumACD_clustered <- coeftest(all_sumACD, vcov = vcovCL, cluster = ~managing_authority)
  all_sumACD_clusterSE <- all_sumACD_clustered[,2]
  all_sumACDE <- lm(GL_pct ~ sum_impact + cast + drill + as.factor(election_id) + as.factor(managing_authority), allprovinces)
  all_sumACDE_clustered <- coeftest(all_sumACDE, vcov = vcovCL, cluster = ~managing_authority)
  all_sumACDE_clusterSE <- all_sumACDE_clustered[,2]
  
  ## Regressions using all provinces, 4 earthquake bin variables
  
  all_sumCD_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                      cast + drill, bin4)
  all_sumCD_1_clustered <- coeftest(all_sumCD_1, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_1_clusterSE <- all_sumCD_1_clustered[,2]
  all_sumCDE_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                       cast + drill + as.factor(election_id), bin4)
  all_sumCDE_1_clustered <- coeftest(all_sumCDE_1, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCDE_1_clusterSE <- all_sumCDE_1_clustered[,2]
  all_sumACD_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                       cast + drill + as.factor(managing_authority), bin4)
  all_sumACD_1_clustered <- coeftest(all_sumACD_1, vcov = vcovCL, cluster = ~managing_authority)
  all_sumACD_1_clusterSE <- all_sumACD_1_clustered[,2]
  all_sumACDE_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                        cast + drill + as.factor(election_id) + as.factor(managing_authority), bin4)
  all_sumACDE_1_clustered <- coeftest(all_sumACDE_1, vcov = vcovCL, cluster = ~managing_authority)
  all_sumACDE_1_clusterSE <- all_sumACDE_1_clustered[,2]
  
  ### Subset data by province to run province-specific regressions
  
  groningen <- allprovinces[allprovinces$contest_name=="Groningen", ]
  drenthe <- allprovinces[allprovinces$contest_name=="Assen", ]
  friesland <- allprovinces[allprovinces$contest_name=="Leeuwarden", ]
  
  groningen4 <- bin4[allprovinces$contest_name=="Groningen", ]
  drenthe4 <- bin4[allprovinces$contest_name=="Assen", ]
  friesland4 <- bin4[allprovinces$contest_name=="Leeuwarden", ]
  
  ## Regressions using Groningen, main earthquake variable
  
  groningen_sumCD <- lm(GL_pct ~ sum_impact + cast + drill, groningen)
  groningen_sumCD_clustered <- coeftest(groningen_sumCD, vcov = vcovCL, cluster = ~managing_authority)
  groningen_sumCD_clusterSE <- groningen_sumCD_clustered[,2]
  
  ## Regressions using Groningen, 4 earthquake bin variables
  
  g_sumCD_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                    cast + drill, groningen4)
  g_sumCD_1_clustered <- coeftest(g_sumCD_1, vcov = vcovCL, cluster = ~managing_authority)
  g_sumCD_1_clusterSE <- g_sumCD_1_clustered[,2]
  
  ## Regressions using Drenthe, main earthquake variable
  
  drenthe_sumCD <- lm(GL_pct ~ sum_impact + cast + drill, drenthe)
  drenthe_sumCD_clustered <- coeftest(drenthe_sumCD, vcov = vcovCL, cluster = ~managing_authority)
  drenthe_sumCD_clusterSE <- drenthe_sumCD_clustered[,2]
  
  ## Regressions using Drenthe, 4 earthquake bin variables
  
  d_sumCD_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                    cast + drill, drenthe4)
  d_sumCD_1_clustered <- coeftest(d_sumCD_1, vcov = vcovCL, cluster = ~managing_authority)
  d_sumCD_1_clusterSE <- d_sumCD_1_clustered[,2]
  
  ## Regressions using Friesland, main earthquake variable
  
  friesland_sumCD <- lm(GL_pct ~ sum_impact + cast + drill, friesland)
  friesland_sumCD_clustered <- coeftest(friesland_sumCD, vcov = vcovCL, cluster = ~managing_authority)
  friesland_sumCD_clusterSE <- friesland_sumCD_clustered[,2]
  
  ## Regressions using Friesland, 4 earthquake bin variable
  
  f_sumCD_1 <- lm(GL_pct ~ sum_impact_1 + sum_impact_2 + sum_impact_3 + sum_impact_4 + 
                    cast + drill, friesland4)
  f_sumCD_1_clustered <- coeftest(f_sumCD_1, vcov = vcovCL, cluster = ~managing_authority)
  f_sumCD_1_clusterSE <- f_sumCD_1_clustered[,2]
  
  ### Placebo test regressions
  
  all_sumCD_p <- lm(GL_pct ~ sum_impact_p + cast + drill, allprovinces)
  all_sumCD_p_clustered <- coeftest(all_sumCD_p, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_p_clusterSE <- all_sumCD_p_clustered[,2]
  all_sumCDE_p <- lm(GL_pct ~ sum_impact_p + cast + drill + as.factor(election_id), allprovinces)
  all_sumCDE_p_clustered <- coeftest(all_sumCDE_p, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCDE_p_clusterSE <- all_sumCDE_p_clustered[,2]
  all_sumACD_p <- lm(GL_pct ~ sum_impact_p + cast + drill + as.factor(managing_authority), allprovinces)
  all_sumACD_p_clustered <- coeftest(all_sumACD_p, vcov = vcovCL, cluster = ~managing_authority)
  all_sumACD_p_clusterSE <- all_sumACD_p_clustered[,2]
  all_sumACDE_p <- lm(GL_pct ~ sum_impact_p + cast + drill + as.factor(election_id) + as.factor(managing_authority), allprovinces)
  all_sumACDE_p_clustered <- coeftest(all_sumACDE_p, vcov = vcovCL, cluster = ~managing_authority)
  all_sumACDE_p_clusterSE <- all_sumACDE_p_clustered[,2]
  
  g_sumCD_p <- lm(GL_pct ~ sum_impact_p + cast + drill, groningen)
  g_sumCD_p_clustered <- coeftest(g_sumCD_p, vcov = vcovCL, cluster = ~managing_authority)
  g_sumCD_p_clusterSE <- g_sumCD_p_clustered[,2]
  d_sumCD_p <- lm(GL_pct ~ sum_impact_p + cast + drill, drenthe)
  d_sumCD_p_clustered <- coeftest(d_sumCD_p, vcov = vcovCL, cluster = ~managing_authority)
  d_sumCD_p_clusterSE <- d_sumCD_p_clustered[,2]
  f_sumCD_p <- lm(GL_pct ~ sum_impact_p + cast + drill, friesland)
  f_sumCD_p_clustered <- coeftest(f_sumCD_p, vcov = vcovCL, cluster = ~managing_authority)
  f_sumCD_p_clusterSE <- f_sumCD_p_clustered[,2]
  
  ### Run main regressions for other parties
  
  # CDA (center right)
  
  ## Create a variable that is vote share for Christian Democratic Appeal, a major center right party
  allprovinces$CDA_pct <- allprovinces$CDA/allprovinces$total_counted
  groningen$CDA_pct <- groningen$CDA/groningen$total_counted
  drenthe$CDA_pct <- drenthe$CDA/drenthe$total_counted
  friesland$CDA_pct <- friesland$CDA/friesland$total_counted
  
  ## Track incumbency for CDA
  allprovinces$CDA_inc <- ifelse(allprovinces$election_id == "TK2017", 0, 1)
  groningen$CDA_inc <- ifelse(groningen$election_id == "TK2017", 0, 1)
  drenthe$CDA_inc <- ifelse(drenthe$election_id == "TK2017", 0, 1)
  friesland$CDA_inc <- ifelse(friesland$election_id == "TK2017", 0, 1)
  
  ## Regressions using all provinces, main earthquake variable
  all_sumCD_CDA <- lm(CDA_pct ~ sum_impact + cast + drill + CDA_inc, allprovinces)
  names(all_sumCD_CDA$coefficients)[names(all_sumCD_CDA$coefficients) == "CDA_inc"] <- "incumbency"
  all_sumCD_CDA_clustered <- coeftest(all_sumCD_CDA, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_CDA_clusterSE <- all_sumCD_CDA_clustered[,2]
  
  ## Regressions using Groningen, main earthquake variable
  groningen_sumCD_CDA <- lm(CDA_pct ~ sum_impact + cast + drill + CDA_inc, groningen)
  names(groningen_sumCD_CDA$coefficients)[names(groningen_sumCD_CDA$coefficients) == "CDA_inc"] <- "incumbency"
  groningen_sumCD_CDA_clustered <- coeftest(groningen_sumCD_CDA, vcov = vcovCL, cluster = ~managing_authority)
  groningen_sumCD_CDA_clusterSE <- groningen_sumCD_CDA_clustered[,2]
  
  ## Regressions using Drenthe, main earthquake variable
  drenthe_sumCD_CDA <- lm(CDA_pct ~ sum_impact + cast + drill + CDA_inc, drenthe)
  names(drenthe_sumCD_CDA$coefficients)[names(drenthe_sumCD_CDA$coefficients) == "CDA_inc"] <- "incumbency"
  drenthe_sumCD_CDA_clustered <- coeftest(drenthe_sumCD_CDA, vcov = vcovCL, cluster = ~managing_authority)
  drenthe_sumCD_CDA_clusterSE <- drenthe_sumCD_CDA_clustered[,2]
  
  ## Regressions using Friesland, main earthquake variable
  friesland_sumCD_CDA <- lm(CDA_pct ~ sum_impact + cast + drill + CDA_inc, friesland) 
  names(friesland_sumCD_CDA$coefficients)[names(friesland_sumCD_CDA$coefficients) == "CDA_inc"] <- "incumbency"
  friesland_sumCD_CDA_clustered <- coeftest(friesland_sumCD_CDA, vcov = vcovCL, cluster = ~managing_authority)
  friesland_sumCD_CDA_clusterSE <- friesland_sumCD_CDA_clustered[,2]
  
  # PVDA (center left)
  
  ## Create a variable that is vote share for Labour Party, the major center left party
  allprovinces$PVDA_pct <- allprovinces$PVDA/allprovinces$total_counted
  groningen$PVDA_pct <- groningen$PVDA/groningen$total_counted
  drenthe$PVDA_pct <- drenthe$PVDA/drenthe$total_counted
  friesland$PVDA_pct <- friesland$PVDA/friesland$total_counted
  
  ## Track incumbency for PvdA
  allprovinces$PVDA_inc <- ifelse(allprovinces$election_id == "TK2012", 0, 1)
  groningen$PVDA_inc <- ifelse(groningen$election_id == "TK2012", 0, 1)
  drenthe$PVDA_inc <- ifelse(drenthe$election_id == "TK2012", 0, 1)
  friesland$PVDA_inc <- ifelse(friesland$election_id == "TK2012", 0, 1)
  
  ## Regressions using all provinces, main earthquake variable
  all_sumCD_PVDA <- lm(PVDA_pct ~ sum_impact + cast + drill + PVDA_inc, allprovinces) 
  names(all_sumCD_PVDA$coefficients)[names(all_sumCD_PVDA$coefficients) == "PVDA_inc"] <- "incumbency"
  all_sumCD_PVDA_clustered <- coeftest(all_sumCD_PVDA, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_PVDA_clusterSE <- all_sumCD_PVDA_clustered[,2]
  
  ## Regressions using Groningen, main earthquake variable
  groningen_sumCD_PVDA <- lm(PVDA_pct ~ sum_impact + cast + drill + PVDA_inc, groningen)
  names(groningen_sumCD_PVDA$coefficients)[names(groningen_sumCD_PVDA$coefficients) == "PVDA_inc"] <- "incumbency"
  groningen_sumCD_PVDA_clustered <- coeftest(groningen_sumCD_PVDA, vcov = vcovCL, cluster = ~managing_authority)
  groningen_sumCD_PVDA_clusterSE <- groningen_sumCD_PVDA_clustered[,2]
  
  ## Regressions using Drenthe, main earthquake variable
  drenthe_sumCD_PVDA <- lm(PVDA_pct ~ sum_impact + cast + drill + PVDA_inc, drenthe)
  names(drenthe_sumCD_PVDA$coefficients)[names(drenthe_sumCD_PVDA$coefficients) == "PVDA_inc"] <- "incumbency"
  drenthe_sumCD_PVDA_clustered <- coeftest(drenthe_sumCD_PVDA, vcov = vcovCL, cluster = ~managing_authority)
  drenthe_sumCD_PVDA_clusterSE <- drenthe_sumCD_PVDA_clustered[,2]
  
  ## Regressions using Friesland, main earthquake variable
  friesland_sumCD_PVDA <- lm(PVDA_pct ~ sum_impact + cast + drill + PVDA_inc, friesland)
  names(friesland_sumCD_PVDA$coefficients)[names(friesland_sumCD_PVDA$coefficients) == "PVDA_inc"] <- "incumbency"
  friesland_sumCD_PVDA_clustered <- coeftest(friesland_sumCD_PVDA, vcov = vcovCL, cluster = ~managing_authority)
  friesland_sumCD_PVDA_clusterSE <- friesland_sumCD_PVDA_clustered[,2]
  
  # SP (left-wing)
  
  ## Create a variable that is vote share for Socialist Party, the major left-wing party
  allprovinces$SP_pct <- allprovinces$SP/allprovinces$total_counted
  groningen$SP_pct <- groningen$SP/groningen$total_counted
  drenthe$SP_pct <- drenthe$SP/drenthe$total_counted
  friesland$SP_pct <- friesland$SP/friesland$total_counted
  
  ## SP were never incumbent
  
  ## Regressions using all provinces, main earthquake variable
  all_sumCD_SP <- lm(SP_pct ~ sum_impact + cast + drill, allprovinces)
  all_sumCD_SP_clustered <- coeftest(all_sumCD_SP, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_SP_clusterSE <- all_sumCD_SP_clustered[,2]
  
  ## Regressions using Groningen, main earthquake variable
  groningen_sumCD_SP <- lm(SP_pct ~ sum_impact + cast + drill, groningen)
  groningen_sumCD_SP_clustered <- coeftest(groningen_sumCD_SP, vcov = vcovCL, cluster = ~managing_authority)
  groningen_sumCD_SP_clusterSE <- groningen_sumCD_SP_clustered[,2]
  
  ## Regressions using Drenthe, main earthquake variable
  drenthe_sumCD_SP <- lm(SP_pct ~ sum_impact + cast + drill, drenthe)
  drenthe_sumCD_SP_clustered <- coeftest(drenthe_sumCD_SP, vcov = vcovCL, cluster = ~managing_authority)
  drenthe_sumCD_SP_clusterSE <- drenthe_sumCD_SP_clustered[,2]
  
  ## Regressions using Friesland, main earthquake variable
  friesland_sumCD_SP <- lm(SP_pct ~ sum_impact + cast + drill, friesland)
  friesland_sumCD_SP_clustered <- coeftest(friesland_sumCD_SP, vcov = vcovCL, cluster = ~managing_authority)
  friesland_sumCD_SP_clusterSE <- friesland_sumCD_SP_clustered[,2]
  
  # VVD (center right)
  
  ## Create a variable that is vote share for People's Party for Freedom and Democracy, a major center right party
  allprovinces$VVD_pct <- allprovinces$VVD/allprovinces$total_counted
  groningen$VVD_pct <- groningen$VVD/groningen$total_counted
  drenthe$VVD_pct <- drenthe$VVD/drenthe$total_counted
  friesland$VVD_pct <- friesland$VVD/friesland$total_counted
  
  ## Track incumbency for VVD
  allprovinces$VVD_inc <- ifelse(allprovinces$election_id == "TK2010", 0, 1)
  groningen$VVD_inc <- ifelse(groningen$election_id == "TK2010", 0, 1)
  drenthe$VVD_inc <- ifelse(drenthe$election_id == "TK2010", 0, 1)
  friesland$VVD_inc <- ifelse(friesland$election_id == "TK2010", 0, 1)
  
  ## Regressions using all provinces, main earthquake variable
  all_sumCD_VVD <- lm(VVD_pct ~ sum_impact + cast + drill + VVD_inc, allprovinces)
  names(all_sumCD_VVD$coefficients)[names(all_sumCD_VVD$coefficients) == "VVD_inc"] <- "incumbency"
  all_sumCD_VVD_clustered <- coeftest(all_sumCD_VVD, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_VVD_clusterSE <- all_sumCD_VVD_clustered[,2]
  
  ## Regressions using Groningen, main earthquake variable
  groningen_sumCD_VVD <- lm(VVD_pct ~ sum_impact + cast + drill + VVD_inc, groningen)
  names(groningen_sumCD_VVD$coefficients)[names(groningen_sumCD_VVD$coefficients) == "VVD_inc"] <- "incumbency"
  groningen_sumCD_VVD_clustered <- coeftest(groningen_sumCD_VVD, vcov = vcovCL, cluster = ~managing_authority)
  groningen_sumCD_VVD_clusterSE <- groningen_sumCD_VVD_clustered[,2]
  
  ## Regressions using Drenthe, main earthquake variable
  drenthe_sumCD_VVD <- lm(VVD_pct ~ sum_impact + cast + drill + VVD_inc, drenthe)
  names(drenthe_sumCD_VVD$coefficients)[names(drenthe_sumCD_VVD$coefficients) == "VVD_inc"] <- "incumbency"
  drenthe_sumCD_VVD_clustered <- coeftest(drenthe_sumCD_VVD, vcov = vcovCL, cluster = ~managing_authority)
  drenthe_sumCD_VVD_clusterSE <- drenthe_sumCD_VVD_clustered[,2]
  
  ## Regressions using Friesland, main earthquake variable
  friesland_sumCD_VVD <- lm(VVD_pct ~ sum_impact + cast + drill + VVD_inc, friesland)
  names(friesland_sumCD_VVD$coefficients)[names(friesland_sumCD_VVD$coefficients) == "VVD_inc"] <- "incumbency"
  friesland_sumCD_VVD_clustered <- coeftest(friesland_sumCD_VVD, vcov = vcovCL, cluster = ~managing_authority)
  friesland_sumCD_VVD_clusterSE <- friesland_sumCD_VVD_clustered[,2]
  
  # PVV (right-wing populist)
  
  ## Create a variable that is vote share for Party for Freedom, the rising right-wing populist party
  allprovinces$PVV_pct <- allprovinces$PVV/allprovinces$total_counted
  groningen$PVV_pct <- groningen$PVV/groningen$total_counted
  drenthe$PVV_pct <- drenthe$PVV/drenthe$total_counted
  friesland$PVV_pct <- friesland$PVV/friesland$total_counted
  
  ## PVV were never incumbent
  
  ## Regressions using all provinces, main earthquake variable
  all_sumCD_PVV <- lm(PVV_pct ~ sum_impact + cast + drill, allprovinces)
  all_sumCD_PVV_clustered <- coeftest(all_sumCD_PVV, vcov = vcovCL, cluster = ~managing_authority)
  all_sumCD_PVV_clusterSE <- all_sumCD_PVV_clustered[,2]
  
  ## Regressions using Groningen, main earthquake variable
  groningen_sumCD_PVV <- lm(PVV_pct ~ sum_impact + cast + drill, groningen)
  groningen_sumCD_PVV_clustered <- coeftest(groningen_sumCD_PVV, vcov = vcovCL, cluster = ~managing_authority)
  groningen_sumCD_PVV_clusterSE <- groningen_sumCD_PVV_clustered[,2]
  
  ## Regressions using Drenthe, main earthquake variable
  drenthe_sumCD_PVV <- lm(PVV_pct ~ sum_impact + cast + drill, drenthe)
  drenthe_sumCD_PVV_clustered <- coeftest(drenthe_sumCD_PVV, vcov = vcovCL, cluster = ~managing_authority)
  drenthe_sumCD_PVV_clusterSE <- drenthe_sumCD_PVV_clustered[,2]
  
  ## Regressions using Friesland, main earthquake variable
  friesland_sumCD_PVV <- lm(PVV_pct ~ sum_impact + cast + drill, friesland)
  friesland_sumCD_PVV_clustered <- coeftest(friesland_sumCD_PVV, vcov = vcovCL, cluster = ~managing_authority)
  friesland_sumCD_PVV_clusterSE <- friesland_sumCD_PVV_clustered[,2]
  
  ### Non-linear model
  
  allprovinces$sum_impact2 <- (allprovinces$sum_impact)^2
  
  quad_CD <- lm(GL_pct ~ sum_impact + sum_impact2 + cast + drill, allprovinces)
  quad_CD_clustered <- coeftest(quad_CD, vcov = vcovCL, cluster = ~managing_authority)
  quad_CD_clusterSE <- quad_CD_clustered[,2]
  quad_CDE <- lm(GL_pct ~ sum_impact + sum_impact2 + cast + drill + as.factor(election_id), allprovinces)
  quad_CDE_clustered <- coeftest(quad_CDE, vcov = vcovCL, cluster = ~managing_authority)
  quad_CDE_clusterSE <- quad_CDE_clustered[,2]
  quad_ACD <- lm(GL_pct ~ sum_impact + sum_impact2 + cast + drill + as.factor(managing_authority), allprovinces)
  quad_ACD_clustered <- coeftest(quad_ACD, vcov = vcovCL, cluster = ~managing_authority)
  quad_ACD_clusterSE <- quad_ACD_clustered[,2]
  quad_ACDE <- lm(GL_pct ~ sum_impact + sum_impact2 + cast + drill + as.factor(election_id) + as.factor(managing_authority), allprovinces)
  quad_ACDE_clustered <- coeftest(quad_ACDE, vcov = vcovCL, cluster = ~managing_authority)
  quad_ACDE_clusterSE <- quad_ACDE_clustered[,2]